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ABSTRACT 

We use large cosmological Smoothed-Particle- Hydrodynamics simulations to study the 
formation and evolution of sub- millimetre galaxies (SMGs). In our previous work, we 
studied the statistical properties of ultra-violet selected star-forming galaxies at high 
redshifts. We populate the same cosmological simulations with SMGs by calculating 
the reprocess of stellar light by dust grains into far-infrared to millimetre wavebands 
in a self-consistent manner. We generate light-cone outputs to compare directly the 
statistical properties of the simulated SMGs with available observations. Our model 
reproduces the submm source number counts and the clustering amplitude. We show 
that bright SMGs with flux > 1 mJy reside in halos with mass of ^ 10^"^ Mq and 
have stellar masses greater than 10^^ M0. The angular cross-correlation between the 
SMGs and Lyman-a emitters is significantly weaker than that between the SMGs and 
Lyman-break galaxies. The cross-correlation is also weaker than the auto-correlation 
of the SMGs. The redshift distribution of the SMGs shows a broad peak at z 2, 
where Bright SMGs contribute significantly to the global cosmic star formation rate 
density. Our model predicts that there are hundreds of SMGs with > 0.1 mJy at 
z > 5 per 1 square degree field. Such SMGs can be detected by ALMA. 
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1 INTRODUCTION 

An array of recent millimetre and sub-millimetre (submm) 
observations revealed the physical properties of sub- 
millimetre bright galajcies (SMG s) with luminosities 
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Tamura et al.' '2009; 'Knu dsen et all 1201 Ol :l Riechers et al.l 
2010 : Hatsukadc ct al._ .201li r" The power sources of such 
luminous SMGs are thought to be very active star forma- 
tion with 100 — lOOOM0yr~^ and/or active galactic nuclei 
(AGNs). Ultra-violet photons from massive stars or AGNs 
are converted to photons in submm to infr ared bands b y dust 
grains via their thermal emission. Alexander et al.l 
find that, in the Chandra Deep Field North survey, X-rays 
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are detected from more than one-third of bright SMGs, 
which indicates a significant contribution from AGNs. 

SMGs provide important information on the star forma- 
tion rate, the stellar initial mass function, and the chemical 
evolution of the galaxies at high redshifts. Because SMGs 
are highly obscured by dust, they often appear dark in UV 
and optical wavebands. Thus measuring their redshifts by 
conventional techniques is very difficult. Moreover, identify- 
ing the optical counterpart is not easy because the spatial 
resolution of the currently operating submm telescopes is 
much worse than that of large optical telescopes. 

So far, most of the SMGs wit h measured redsh i fts ar e 
at z < 3 (|Swinbank et al.l |2004| : IChapman et all 120051 ) . 
but a few SMGs have been found a t higher redshifts 
(ICapak et al.l l2008l: ICoppin et all l2009l : iDaddi et al.] l2009l : 
iKnudsen et al.l 2010l). SMGs generally show strong cluster- 
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2 Shimizu, Yoshida and Okamoto 



the 350 fim and 500 selected SMGs in the Herschel- 
ATLAS survey. The SMGs are strongly clustered while the 
250 selected sample showed weak or no clustering sig- 
nals. It is thought that SMGs are very massive systems 
with large gas reservoir 10^^ Mp) and with large stel- 
lar masses lO" Mr) (|Greve et al.l l2005l : iTacconi et"al] 
l2006l : IWardlow et al.llioill '). Also, SMGs are likely to be an- 
cestors of m assive elliptical galaxies in the local universe 
l|Lillv et al.llT996: Small et al. 2004 ) . 

IChapman erahl (|2005l ') report that the redshift distri- 
bution of SMGs shows a peak at z ~ 2, similarly to that 
of AGNs, although the size of the sample is small. Several 
SMGs have counterparts of star-forming galaxies such as Ly- 
man break galaxies (LBG s ) and Lyman a em i tters (LAEs) 
jGeach et al.l |2005| . l2007l : iSeelen et al.] l2008l : iDaddi et al.l 
I2OO9I ). However, many of high-redshift galaxies are not seen 
as SMGs probably because their submm fluxes are below 
typical observational flux limits (> 1 mjy). Because ob- 
servations so far detect very bright SMGs, it remains un- 
clear if SMGs are ge nera lly associated with LBGs/LAEs. 
iMatsuda etHI l|20oi ) and lTamura et all (|2010l ') failed to de- 
tect submm sources at the location of LAEs in a protoclus- 
ter region even though the target galaxies are very bright in 
Lya- 

iDaval et al.l (|2010l l use cosmological simulations and 
found that the submm fluxes of high-z LAEs would 
be less than 0.1 mJy. They argue, however, that 
many of such galaxies c an b e detected by ALMA. 
lYajima. Umemura. fc Moril ()201lh study the evolution of 
submm flux of a simulated galaxy by using detailed three- 
dimensional radiative transfer calculations. They show that 
the simulated galaxy in a very early LAE phase can be de- 
tected by ALMA. Observing galajcies in the early formation 
phase is important to understand the early chemical evolu- 
tion. 

One of the most important observational quantities 



is the source number count of SMGs( 
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20091. I2OIOI: IScott et all I2OIOI: iHatsukade et al.ll201ll). Un- 



fortunately, ther e remain sub s tantia l uncertainties in the- 
oretical models. iBaugh et al.l l|2005h calculate the source 
number count of SMGs using a semi-analytic galaxy for- 
mation model. Interestingly, in order to reproduce the ob- 
served source number count, they make a crucial assump- 
tion that the stellar initial mass fun ction (IMF) is 'to p- 
heavy' when s tarburst occurs ( see al so lLacev et al.l l|2010l) ). 
More recently, iFontanot et al.l l|2007h reproduced the source 
number count of SMGs without assuming to p-heavy IMF. 
They argue that the major difference between iBaugh et al.l 
1120051) and their stud y may be in the gas cooling model. In 
iFontanot et al.l (|2007l ). the h ot gas in a galaxy is assumed to 
cool more efficiently than in lBaugh et al.l ll2005h. However, 



becau se of the efficient cooling, the model of Fontanot et al.l 
l|2007l ) overproduces massive galaxies at lower redshifts. It 
appears that some efficient star-formation recipe is needed 
in the semi-analytic models, in order to reproduce the source 
number count of SMGs. It is interesting to see if the same 
is true for cosmologica l hydro dynamic simulations of galaxy 
formation. iDave et al.l l|2010l ) study the physical properties 
of SMGs using cosmological simulations. They claim that 
smooth infalling gas or/and gas-rich satellites are impor- 
tant to produce large star formation rates. Although their 



finding seems reasonable, they employ a simple model that 
galaxies with very large star-formation rates are identified 
as SMGs. They do not calculate dust absorption and the re- 
sulting submm flux. It is important and timely to study in 
detail the statistical properties of SMGs over a wide range 
of redshift using cosmological hydrodynamic simulations by 
considering dust absorption and re-emission consistently. 

In this paper, we study the statistical properties of 
SMGs such as the stellar mass function, the source number 
count and the angular correlation function. To this end, we 
perform large cosmological hydrodynamic simulations based 
on the standard ACDM cosmology. Our simulations follow 
star formation, supernova feedback, and metal enrichment 
self-consistently. For the galaxies identifled in our cosmo- 
logical simulations, we calculate the spectral evolution and 
dust extinction to esti mate the FIR luminosit y. In our ear- 
lier work (Shim izu, Yo shida fc Okamotci|l201ll ). we used the 
same set of simulations to study the statistical properties of 
LAEs at z = 3.1. Our model successfully reproduced all the 
available observational data at 2; = 3.1 such as the Lyman-Q 
luminosity function, the angular correlation functions, and 
the Lyman-Q equivalent width distribution. Our aim in the 
present paper is to build a consistent theoretical model of 
galaxy formation that reproduces the observed properties of 
both ultra-violet selected galaxies and SMGs. 

Throughout the present paper, we adopt the ACDM 
cosmology with the matter density f^M = 0.27, the cosmo- 
logical constant = 0.73, the Hubble constant h — 0.7 in 
units of Ho = 100 km s~^ Mpc~^ and the baryon density 
Qb = 0.046. The matt er density fluctuatio ns are normalised 
by setting erg = 0.81 (ISpergel et al.ll2003l ). All magnitudes 
are expressed in the AB system, and the Lya EWlyq values 
in this paper are in the rest frame. 



2 THEORETICAL MODEL 

Our simulation code is based on an updated version of 
the Tree-PM smoothed particle hydrodynamics (SPH) code 
GADGET-3 whlch is a s uccessor of Tree-PM SPH code 
GADGET-2 ISpringelllioOsI ). We employ iV = 2 x 640^ par- 
ticles in a comoving volume of 100 Mpc on a side. The 
mass of a dark matter particle is 2.41 x 10*/i~^Mq and 
that of a gas particle is 4.95 x 10^/i~^Mq, respectively. We 
implement relevant physical processes such as star forma- 
tio n, supernova feedback and chemica l enrichment follow- 
ing lOkam oto, Nemmen fc Bowei' ('2008') , 'Okamoto, & Fren^ 
pOOal and Okamoto ot al. (2010). In Okamoto ct al. (2010i), 
two wind models for supernova (SN) feedback are ex- 
amined. One model assumes that the initial wind speed 
of gas around stellar particles that trigger SN is pro- 
portional to the local velocity dispersion of dark mat- 
ter (o;)^_The_model is motivated by recent observations 
of iMartinI l|2005l ). The other model assumes a constant 
initial wind speed, as has be en used in many cosmologi- 
cal hydrodynamic simulations |Springel fc Hernauistll2003al : 



iNagamine. Springel fc Hernguistl |2004| ). The major differ- 
ence between the two models is the wind mass-loading per 
star formation rate. In the former model, the mass-loading 
is proportional to w hereas it is cons t ant fo r all galax- 
ies in the latter model. lOkamoto etHI l|2010l ) show that 
the former model reproduce the observed satellite luminos- 
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ity function and the luminosity-metallicity relation in the 
Local Group satellites. We set the initial wind speed of 
gas to be five t i mes t he local velocity dispersion, following 
lOkamoto et"all (|201(]| ). 

In our previous paper l|Shimizu. Yoshida fc Okamotd 

we studied the statistical properties of UV-selected 
galaxies. We implemented wavelength dependent dust ab- 
sorption to calculate the UV and Lyman-Q luminosities of 
the galaxies. In the present paper, in order to calculate the 
FIR luminosity, we assume that the absorbed UV photons 
are re-emitted in FIR by thermal dust grains. 



2.1 Dust absorption 

The treatment of dust absorpt ion is described in 
IShimizu. Yoshida fc Okamotd |20l3). Briefly, we calculate 
the optical depth Td(A) for UV continuum photons on the 
assumption that dust is distributed in a sphere with a cer- 
tain radius around a galaxy. Then we use only one length 
scale, the effective radius of dust distribution r^, to evaluate 
the optical depth. 



3Sdg(A) 
4adS ' 



(1) 



where and s are the dust grain size and the material den- 
sity of dust grains, respectively, and Ed = Md/{Tvr^) is the 
dust surface mass density. For simplicity, we assume a single- 
sized dust component. We use the d ust optical constant Q {\) 
as a function of wavelength given in lDraine fc Led(|l984l ). so 
that we obtain realistic SEDs for individual galaxies. We set 
Od = 0.1 fim and s = 2.5 g cm~'^ f or dust grain s that are pro- 
duce d by supernovae l|Todini fc Fe rarra 200ll : iNozawa et al.l 
I2OO3I ). In this case, the resulting extinction curve is flatter 
than the SMC extinction curve and the Calzetti extinction 
curve, because Q(A) for large grains is ne arly independent 
of wavelength (|Hirashita et al ] |2008l . l20ld ). Intriguingly, re- 
cent observations of ultraluminous infrared galaxies at z ~ 1 
suggest that the extinction curve is fl atter than the Calzetti 
extinction curve (|Shimizu et al.ll201ll '). Furthermore, the ex- 
tinction curve for the SN-prod uced dust size distribution re - 
produces the observation well (|Hirashita et al.ll2008l . [2OI0I ). 

Although our single-sized dust model might seem crude, 
the resulting extinction curve does not significantly affect 
our main results presented in the paper. In section 4, we ex- 
plicitly show that adopting the steep SMC extinction curve 
yields virtually the same statistical quantities for the SMG 
population. 

We calculate the dust surface mass density as Ed = 
Afd/(7rrd) where Afd and rd are the dust mass and the ef- 
fective dust distribution radius of each simulated galaxy, re- 
spectively. The dust mass is calculated by 



Md = O.Olmj, — 



(2) 



where mg, Z and Zq denote the gas mas s, the gas metallicit y 
and the solar metallicity, respectively (jPraine et al.|[2007l ). 
We adopt Zq — 0.02, and then the above equation gives 
Md — O.Qlrrig for Z = Zq, which means that one percent 
(in mass) of the gas in a simulated galaxy with Z — Zq is 
locked in dust grains. Hence the dust mass scales with the 
gas metallicity. We assume that the effective dust distribu- 
tion radius of a galaxy rd is a fraction of its virial radius. 
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Figure 1. UV (1500 A at rest-frame) luminosity function of the 
simulated galaxies at 2 = 2.5. The solid line shows our simulation 
result, which we compare with th e observational data (points with 
error bars) of lOesch et"all 1I2OIOI') . 



rd = e X rvir 



(3) 



where e is a global parameter to be kept constant for all the 
galaxies. The numerical value of e is determined by matching 
the UV luminosity function at z = 2.5 (see Fig.[l}. The best 
fit value at z — 2.5 is e = 0.09. 

We calculate the escape fraction of the UV continuum 
photons fc using the following equation: 

1 _ p-^d(A) 



/c(A) = 



^d(A) 



(4) 



The FIR luminosity Lfir is then set equal to the total energy 
of UV photons absorbed by dust. 



(5) 



where L'"'(i^) is the intrinsic luminosity (per frequency) of 
the galaxy, and UJ'^^ is the luminosity after dust absorp- 
tion. The above integration is evaluated over from UV to 
optical wavelengths. In the next subsection, we describe the 
calculation of the FIR fiux. 



2.2 Calculation of the FIR Flux 

Wo calculate the FIR flux following iHirashita fc Ferraxd 
(2002). For the single-sized dust model, the monochromatic 
luminosity at frequency u is written as 



L^"=* = 47rMdK.B.(rd) 



(6) 



where k^, Bv and Td are the absorption coefficient, the 
Planck function and the dust grain temperature, respec- 
tively. The absorption coefficient Ku in FIR is well described 
by a power-law — ko{v/uo)^ with /3 = 1 ~ 2. 
Then the total FIR luminosity is 



Ldust n 
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Figure 2. The spatial distribution of our simulated SMGs on the past light cone of an observer. We plot SMGs with submm flux greater 
than 0.1 mjy at 270 GHz. The redshift range is from z = 1.3 to z = 8.5. The field of view of the light cone is 40 arcmin on a side. The 
colour indicates the submm flux; blue for faint SMGs and yellow to red for bright SMGs. Also the point size scales with the submm flux. 
The smallest points and the largest points correspond to 0.1 mJy and for 30 mJy, respectively. 



fcTd 



-dx. 



(7) 



where k, hp and c are the Boltzmann constant 



the Planck 

constant, and the speed of light, respectively. We denote 
X = hpV / (kTg^). We set a normalization parameter vq = 
3.0 X lO'^^Hz which corresponds to the wavelength of lOO^m. 
The absorption coefficient kq is given by kq ~ SQi/q (4as)~^, 



where Q^^ is the dust optical constant at vo- The value of 
kq is 46.95 g~^cm^. Note that Q^^/a, and hence kq, are 
independent of dust size. 

Using equation (0, we derive the dust temperature Td 

as 



Td = 7.64 



Md/Mc 



1/6 



(8) 
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where we set /3 = 2 assuming graphite and silicate dust 
grains IjPraine fc Tei[l984). Finally, the observed submm 
flux of a galaxy is obtained from 

robs _ (1 + + , , 

" 47vdl ' 

where di^ is the luminosity distance to the galaxy. The dust 
temperature for our sample galaxies is typically Td — 33 K. 
To directly compare our model with the observational re- 
sults of the AzTEC observation, we consider the filter re- 
sponse function of the AzTEC detector when we calculate 
the submm flux at the observed frame 1.1 mm. 



3 RESULT 

3.1 UV luminosity function 

We flrst examine whether our model reproduces the UV lu- 
minosity function. This is an important check because the 
total FIR luminosity is equal to, by definition in our model, 
the total UV luminosity absorbed by dust grains. Fig. [1] 
shows the UV (1500 A at rest-frame) luminosity function 
of our model at z = 2.5. Clearly, our model reproduces 
the observational data very well. Given the equally good 
agreement between our model predi ctions for z = 3.1 LAEs 
ish imizu. Yoshida fc Okamotol201lh and other available ob- 
servational data of UV-selected galaxies (Shimizu et al., in 
preparation), we are confident that our cosmological sim- 
ulations can be used to build a self-consistent model for 
UV-selected galaxies and SMGs at high redshifts. 
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Figure 3. The source number counts of SMGs for six observed 
frame fluxes. The thick solid line shows the number count at 
the observed frame 1.1 mm (AzTEC detector band). The thick 
dash-dotted line is for the observed frame 850 fim, and the other 
lines are for the observed frames 3.0 mm, 1.4 mm, 730 /im and 
450 fira which correspond the five ALMA bands, respectively. The 
filled circles with error bars show the obs e rved number counts 
at 1.1 mm b y AzTEC ( Greve et al l '2004"; 'Laurent et"al] |2005|: 
Perera et aP I2OO 8: Aust crmann et al. 2009, 2010: IScott et al.l 
2OIOI : iHatsukade et al.ll2oTlli . The filled triangles with error bars 



show the observed number counts at 850 fim (iChapman et all 
I2OO5I) . 



3.2 The light cone output 

Most of the recent observations do not provide redshifts for 
SMGs. In order to directly compare our model with the ob- 
servations, we generate a light cone output which extends 
from z = 1.3 to z = 8.5 using a number of simulation out- 
puts. We coordinate the output redshifts to fill the volume 
of a lightcone from z = 1.3 to 2; = 8.5 without gaps. We 
then randomly shift each simulation box so that the same 
objects (at different epochs) do not appear multiple times 
on a single line-of-sight. Fig.[2]shows the distribution of sim- 
ulated SMGs on the light cone from z — 1.3 to 2: = 8.7. For 
this figure, we set the SMC detection limit to be 0.1 mJy 
at 270 GHz. Most of the SMGs are located at z = 1 - 4, 
but there are dozens of SMGs even at z > 6. Such faint 
SMGs can be detected by ALMA with less than one hour 
integration using 16 antennas at the 350 GHz band. The 
high redshift SMGs will provide invaluable information on 
star formation activities in the early universe. 



3.3 The source number counts 

We use the generated light cone output to conduct mock ob- 
servations of the SMGs. Fig. [3] shows the number counts of 
SMGs. The filled circles and the filled triangles are the obser- 
vational results at the observed frame 1.1 mm and 850 fim, 
respectively. The lines shows our model predictions for var- 
ious bands. To directly compare our model with the obser- 
vation, we use the filter response function of the AzTEC 



detector when we calculate the submm fiux at the observed 
frame 1.1 mm. 

Our model reproduces the observational results remark- 
ably well. The 850 /im source count is somewhat smaller 
than the observation at S < 1 mJy, but the error bars are 
also large at the faint end. The thin solid, dashes, dotted and 
double-dotted lines are the predicted SMG number counts 
at observed frames 3.0 mm, 1.4 mm, 730 nm and 450 /im, 
which correspond to the ALMA detector bands. The number 
of SMGs increases with shorter wavelength, suggesting that 
SMGs are detected easier at shorter wavelengths; this is an- 
other manifestation of the negative k-correction for SMGs. 

3.4 The SMG number density evolution 

In this subsection, we study the redshift evolution of the 
SMG population in more detail. Fig. |4] shows the evolution 
of the number density of SMGs. The dashed and dash-dotted 
lines represent the number counts for the detection limit of 
0.1 mJy and 1 mJy at the observed frame 1.1 mm. We also 
plot our model number count at 850 ^m with the detection 
limit 5 mJy. The triple- dotted line show s the o bservational 



result at 850 /jm from Chapman et akl l|2005l ). Note that 



the sample of IChapman et aL I (|2005l ) consists of SMGs in 
multiple fields with various areas. Thus we normalize the 
observed redshift distribution arbitrarily to make the com- 
parison easy. We also plot the number of all the simulated 
galaxies (the solid line) with halo mass more than 10^° Mq. 
Regardless of the selection, the SMG number counts peak 
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Figure 4. The evolution of the number density of SMGs as a 
function of redshift. The dashed and dash-dotted Unes represent 
our model prediction for the detection limit of 0.1 mjy and 1 mjy 
at the observed frame 1.1mm, respectively. The dotted line is also 
our model prediction for a high detection limit of 5 mjy at the 
observed frame 850 fim , to be compared with the observation of 
IChapman et al.l l l2005^ . For reference, we plot the number density 
of all simulated galaxies (solid line) with halo masses greater than 
10^° Mq. 

at z ~ 2 — 3. Our model lacks very bright SMGs at high 
redshifts. This is likely due to the small simulation box size 
(cosmic variance). Note also that, in general, determining 
the redshift distribution of the SMGs are not trivial if only 
submm observations are used. 

In our model, approximately 90% of SMGs with > 
0.1 mJy at 1.1 mm are at 2 > 2. This implies that most 
of the observed SMGs are high-z galaxies, and the contribu- 
tion from low-2 galaxies is small. 

We calculate the contribution from SMGs to the global 
cosmic star formation rate (SFR) evolution. Fig.[5]shows the 
sum of the SFRs of all the simulated galaxies (solid line), 
that of SMGs with > 0.1 mJy at observed frame 1.1 mm 
(dashed line) and that of SMGs with > 1 mJy at observed 
frame 1.1 mm (dot-dashed line). 

The contribution of SMGs to the cosmic SFR is small 
at z > 6. However the contribution rapidly increases with 
decreasing redshift. At 2 --^ 2, Bright SMGs with 5 > 1 mJy 
at 1.1 mm contributes roughly one third to the total SFR. 
Clearly, it is important to include SMGs to estimate the true 
global star-formation rate at 2 ~ 1 — 4. 

3.5 The stellar mass of SMGs 

Fig. |6] shows the SFR (the left panel) and specific SFR (the 
right panel) as a function of the stellar mass. We plot all the 
SMGs in our light cone data from z = 1.3 to 2; = 8.7. There 
are tens of SMGs with the SFR greater than 1000 M0yr"\ 
whereas the specific SFR is roughly constant (~ 10~^ yr"'^)- 
Bright SMGs have very large stellar masses in excess of 
10^^ Mq. Some of our simulated SMGs with largest stel- 
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Figure 5. The star formation rate density contributed by SMGs. 
The solid, dashed and dash-dotted lines are for all the galaxies, 
SMGs with S > 0.1 mjy at observed frame 1.1 mm and SMGs 
with S > 1 mjy at observed frame 1.1 mm, respectively. 



lar masses (> 10^^ Mq) consist of multiple clumps within a 
single host halo. However, we do not treat each such clump 
of star particles as an individual galaxy. Essentially we re- 
gard one halo as a SMG. This crude definition is motivated 
by the fact that the current observations do not have high 
enough angular resolution to resolve the internal structure 
within a halo. See also the discussion in section 3.6 and Fig. 
[To] which shows the host halo mass of the SMGs. 



3.6 The distribution of SMGs and the angular 
correlation functions 

Fig. [7] shows the projected distribution of the SMGs on the 
lightcone. We plot the SMGs with S > 0.1 mJy at the ob- 
served frame 1.1 mm. As is naively expected, brighter SMGs 
are strongly clustered, whereas faint SMGs appear to be dis- 
tributed more or less homogeneously over the field of view. 

Next, we calculate the angular two-point correlation 
function (ACF) of the SMGs to quantify the clustering 
strength. In Fig. [S] the solid line is the correlation function 
of our simulated SMG sample. For comparison, we have cho- 
se n the same detection limit of 1.26 mJy as the observation 
of lWilliams et al.l (|201lf ). Our simulation result is consistent 
with the observational data although with substantial scat- 
ters for both the data. 

We calculate the angular cross correlation between the 
simulated SMGs and other UV-selected galaxies. We use one 
simulation box at 2 = 3.1 for this calculation. We consider 
two fiux limits for our simulated SMGs, 0.1 mJy and 1 mJy 
at observed frame 1.1 mm. The number of simulated SMGs 
for the two cases is then 993 and 18, respectively. We identify 
simulated galaxies which sati sfy EWLyg > 30 A and L°%^„ > 
1.0 X 10*^[erg/s] as LAEs ((Sh imizu. Yoshida fc Okamotd 
I2OIII ). Here ^Ly^ and EWLya are the observed Lya lumi- 
nosity and the rest frame Lya equivalent width. We define 
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Figure 6. The star formation rates of the simulated SMGs (left) and the specific star formation rates (right) as a function of the stellar 
mass. We use bigger size points for brighter SMGs. Also the colour indicates the submm flux. To reduce the number of points, we 
randomly chose one tenth of the galaxies for stellar masses below 10^^ Mq. 
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Figure 7. The projected distribution of SMGs in a field-of-view 
of 40 arcmin X 40 arcmin. The left and right panel show SMGs 
with S > O.l mJy. The point size scales with the submm flux. 
The smallest points and the largest points are for ~ 0.1 mJy and 
for ~ 30 mJy, respectively. Also the point colour indicates the 
submm flux. 

LBGs at 2: ~ 3 by the following selection: 

23 < ZAB s; 26.0, (10) 
Uab - Vab > 0.8, Vab - ZAB < 2.7, (11) 

Uab - Vab > 1.8 {Vab - ^ab) + 1.6, (12) 
where Uab, Vab and zab denote AB magnitude of U band. 



V band and z band, respectively. We select the LBGs in 
this way consistently with the observation of 'Yoshida et ajj 
(2008 ). Finally, we locate 2913 LAEs and 3211 LBGs in the 
simulation box at z ~ 3. 

Fig. M shows the SMG-LAE and SMG-LBG angular 
cross-correlation functions. For comparison, we also plot the 
auto correlation function of SMGs (solid circles with error 
bars) and the cross-correlation between SMGs and dark ha- 
los with mass greater than > 10^" M© (stars with error bars). 
For the SMG flux hmit of 0.1 mJy, the SMG-LBG cross cor- 
relation is significantly stronger than the SMG-LAE cross 
correlation. The former is indeed very similar in the shape 
and the amplitude to the auto-correlation of the SMGs, im- 
plying that the SMGs are a similar population to the LBGs 
(see also the SFR shown in Fig. |6]). We have checked that 
many UV bright galaxies in our sample are identified as 
SMGs. 



3.7 The host halo mass of SMGs 

Fig. 1 101 shows the submm flux of simulated SMGs as a func- 
tion of their host halo mass. We use all the SMGs in our 
light cone data from z = 1.3 to z = 8.7. The submm flux of 
simulated SMGs roughly scales with their host halo mass, 
but with large dispersions for low mass objects. We show the 
dispersion of host halo mass and submm flux in each mass 
bin. We also plot dark matter halo mass function at z — 2. 
Clearly, bright SMGs with {S > 1 mJy) reside in rare, large 
halos with several times 10^^ M©. 

It is important to note that, because the angular resolu- 
tion of the current mm/submm observations is still limited 
to ~ 30 — 60 arcseconds, multiple galaxies in a galaxy-group 
size halo are not well resolved. Future observations using 
ALMA will likely flnd multiple sources in massive halos. 
/.From Fig. 1101 we expect that SMGs hosted by smaller halos 
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Figure 8. We plot the angular two-point correlation function of 
the simulated SMGs (solid line with error bars). The points wit h 
error bars are the observational data of IWilliams et al] l l201l[) . 
Different symbols refer to different catalogues of SMGs. The circle 
and triangle points are S.Sc sources and 3(t sources, respectively. 



with several times 10 Mq can be detected by observation 
with a detection limit of 0.1 mjy. 



4 CONCLUSIONS AND DISCUSSION 

We have studied a variety of statistical properties of 
SMGs using a large cosmological hydrodynamic simu- 
lation. The simulation follows the formation and evo- 
lution of star-formi ng galaxies by emplo ying a new 
feedback model of |Okamoto et al.l (|2010l ). Earlier in 
IShimizu. Yoshida fc Okamotol ( 2011 ). the same simulation 
was used to study the properties of LAEs at z = 3.1. We 
have shown in the present paper that our SMG model is con- 
sistent with a number of observational data currently avail- 
able. The redshift distribution of the SMGs peaks broadly 
at 2 ~ 2. Most of the SMGs are at 2 > 2, and the contribu- 
tion from low-2 (2 < 1) galaxies to the submm source count 
is small. The contribution from the SMGs to the global star 
formation rate density increases with decreasing redshift. 
Bright SMGs contributes significantly to the star formation 
rate at 2 ~ 2 — 3. We calculate the angular two point cor- 
relation to quantify the clustering amplitud e of SMGs. Our 
result is consistent with the observations of lWilliams et al] 
l|201ll) . We also study the angular cross correlation func- 
tions between SMGs and LBGs and the correlation between 
SMGs and LAEs. The former is significantly stronger than 
the latter. Considering the high star-formation rates of the 
simulated SMGs, we argue that SMGs are similar population 
to LBGs. Finally, we explicitly showed that the bright SMGs 
preferentially reside in massive halos (> lO^'^ Mq) and that 
their typical stellar mass are greater than 10^^ M©. 

In our simple model of a single-sized, supernovae pro- 
duced dust, the resulting extinction curve is flatter than the 
standard ones such the SMC extinction curve or Calzetti ex- 
tinction curve. This mean that the reprocess of UV photons 



to FIR by dust works more efficiently in our model. However, 
interestingly, we have found that the characteristic feature 
of the extinction curve does not significantly affect our main 
results. 

For direct comparison, we calculate the FIR luminosi- 
ties for our galaxy samples using the SMC extinction curve. 
We then estimate the SMG number count following oth- 
erwise the same procedures as in our fiducial model. Fig. 
111! compares the source number count of our model (solid 
line) and the one calculated using the SMC extinction curve 
(dashed line) . Clearly the difference between our base model 
and the other is very small. This can be explained by the 
fact that we normalize the overall luminosity of the galax- 
ies by matching the UV luminosity function at rest-frame 
1500 A to the observed one (see Fig. [1} • Even for a differ- 
ent extinction model, the overall level of extinction in UV 
is always normalized at 1500 A. Therefore, the difference in 
the exact shape of the extinction curves near UV does not 
affect significantly the SMGs number count in our model, as 
is explicitly shown in Fig. 1111 Note, however, that the UV 
to optical colour of individual galaxies is visibly affected 
by the extinction curve. Our base model predicts that sim- 
ulated galaxies appear bluer in UV to optical range than 
in the case with SMC or Calzetti extinction curve mod- 
els. Fig. [12] shows, as an illustrative example, the SED of 
one of our SMG samples from UV to millimetre. The SED 
for the same galaxy but with SMC extinction curve is also 
shown. We compare them with the observed S EDs for several 
SMGs (jMichalowski. Hiorth fc WatsonI lioiol ) . Some SMGs 
have flat (blue) SEDs at UV to optical range whereas oth- 
ers appear redder. The full SEDs are measured only for a 
limited number of samples. It would be highly interesting to 
observe SMGs in detail in rest frame optical to UV in order 
to derive the dust extinction law, and hence the dust size 
distribution, for star burst galaxies. 

Next, we discuss the source confusion limit of submm 
observations. We have seen in previous sections (Fig. [6] and 
Fig. llOf) . that some simulated SMGs have large host halo 
masses with > lO^'^M© and/or large stellar masses with 
> IO'^Mq. Such massive systems often have multiple sub- 
structures (main and satellite galaxies). Basically each such 
satellite galaxy should be treated as an independent galaxy. 
Note, however, that the angular resolution of the current 
submm/mm observation is worse than that of optical obser- 
vation. For example, the resolution of the AzTEC telescope 
is 30 arcseconds which corresponds 260 kpc (physical scale) 
at 2 = 2. The virial radii of the galactic halos in our simula- 
tion are only as large as ^ 200 kpc (physical size) at 2 = 2. 
Thus, our treatment that a halo hosts one SMG as a whole 
is reasonable, even realistic, if we compare with observations 
with angular resolution of 0.5-1 arcminutes. 

Finally, we discuss the evolution of the stellar mass 
function. Our galaxy formation model is deflned at high 
redshift, rather than being calibrated against galaxies at 
the present epoch. We choose a few model parameters 
such as the overall normalization of dust extinction to re- 
produce the observational data o f star-form ing galaxies at 
high-2 {z ~ 3) ( Shimizu. Yoshi da fc Okamot o 2011 ). Con- 
trastin gly, previous sem i -analy tic models (e.g. jBaugh et al.l 
(|2005l ) and lLacev et al.l (|2010l )) determined the main phys- 
ical parameters to match t he observatio n s of t he present- 
day galaxies. For example, iBaugh et al.l (|2005l ) needed to 
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Figure 9. We plot the angular cross-correlation with two other populations of galaxies; the SMG-LAE cross-correlation (squares with 
error bars) and the SMG-LBG cross-correlation (triangles with error bars). For comparison, we also plot the angular auto-correlation of 
the SMGs (solid circles with error bars) and the angular cross correlation between the SMGs and dark matter halos with 10^" (stars 
with error bars). We set two detection limits for the submm flux; 0.1 mJy for the left panel and 1 mJy for the right panel, respectively. 



10 I 1 1 1 1 1 1 1 1 




Host iialo mass [M,^^] 



Figure 10. The mean submm flux of our simulated SMG samples 
as a function of their host halo mass. We also plot the halo mass 
function at 2; = 2 (the solid line). Bright SMGs with 5 > 1 mJy 
are hosted by dark halos with mass greater than ~ 10^^ Mq. 
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Figure 11. SMG number count at 1.1 mm for the two dust ex- 
tinction models; our model based on a single-sized dust (solid) 
and the other using the SMC extinction curve (dashed). 



adopt an extreme IMF (top-heavy IMF) for starburst galax- 
ies at high redshift to re produce the observe d SMG number 
counts. More recently. Fon tanot et al.l l|2007h reproduced the 
SMGs number count adopting the Salpeter IMF for all their 
galaxy samples. However, their model a ssumes more efficient 
gas cooling in g a lactic haloes than in iBaugh et all (|2005l ) 
and iLacey et al. (|2010l ) . Consequently, very massive galax- 
ies with high star- formation rates are formed, which results 
in disagreement in the stellar mass function at low redshifts 
{z < 1). Clearly it is important to study the stellar mass 
function of our model. 

Fig. [13] shows the stellar mass function of our simulated 



galaxies at z = 1. We also plot the observational data (solid 
points with error bars) of Mortlock ct al. (2011). Our model 
is in good agreement with the data from the GOODS NIC- 
MOS survey as is clearly seen in Fig. [TS] This provides yet 
another support for our galaxy formation model that covers 
a multiple populations, from LAEs, LBGs to SMGs. 

In order to study the present-day stellar mass function, 
we have run our cosmological simulation down to z = 0. 
Our simulation over-predicts the number density of mas- 
sive galaxies with stellar mass greater than IO'^^Mq. There- 
fore, although we have successfully constructed a consistent 



10 Shimizu, Yoshida and Okamoto 




1000 



10 



10° 10° 
wave length [A] 



10' 



Figure 12. An example SED of one of our simulated SMGs. The 
thick solid line and thin dashed line represent the dust-extinct 
SED for our simple dust model and that for the SMC extinc- 
tion curve, respectively. We also show the observ ed SEDs for sev- 
eral SMGs l|Michalowski. Hiorth fc Watsonll201(t by open circles, 
open squares, open triangles and open diamonds. The luminosity 
of the observed galaxies are arbitrary normalized for this compar- 
ison. 
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Figure 13. The stellar mass function at 2 = 1. The solid 
line is our simulation result. The points with error bars show 
the observationjil re sult from the GOODS NICMOS survey 
llMortlock et al.ll201ll) . 



velocity dispersions are greater than 141 km/sec 0. Then 
the resulting stellar mass function closely agrees with the 
observational data at z = derived by Li & White (2009). 
We argue that the kind of strong feedback, in terms of star 
formation efficiency, may be needed for viable galaxy forma- 
tion models, in order for them to reproduce all the available 
data from low to high redshifts. Constructing such a perfect 
model is the ultimate goal of the study of galaxy formation, 
but it is beyond the scope of the present paper. Further stud- 
ies on the efficiency of star formation in galaxies are clearly 
needed using multiple approaches. 

Together with our previous study on Lyman - q emit- 
ters at z = 3. imizu, Yoshida fc Oka moto 2011), we now 
provide a unified galaxy population model within the stan- 
dard A Cold Dark Matter cosmology, which reproduces si- 
multaneously the statistical properties of UV-selected star- 
forming galaxies and sub-millimetre galaxies at high red- 
shifts. 
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